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Abstract 

The general problem of two-phase transport in phase-field models is 
analyzed: the flux of a conserved quantity is driven by the gradient of 
a potential through a medium that consists of domains of two distinct 
phases which are separated by diffuse interfaces. It is shown that the 
finite thickness of the interfaces induces two effects that are not present 
in the analogous sharp-interface problem: a surface excess current and 
a potential jump at the interfaces. It is shown that both effects can be 
eliminated simultaneously only if the coefficient of proportionality between 
fiux and potential gradient (mobility) is allowed to become a tensor in the 
interfaces. This opens the possibility for precise and efficient simulations 
of transport problems with finite interface thickness. 

1 Introduction 

Phase-field models have recently enjoyed a rapidly growing popularity as a com- 
pact and elegant simulation tool for moving boundary problems in such diverse 
fields as solidification (TJ [5], fluid dynamics [3] or solid-state transformations 
[H [S] . Their technical advantage resides in the implicit representation of inter- 
faces by one or several phase fields, i.e. fields that are defined in the entire space, 
take constant values within the bulk of each domain, and exhibit smooth but 
steep variations in well-localized interfacial regions. The tedious procedure of 
front tracking is avoided by introducing equations of motion for the phase fields 
that are coupled to the relevant transport variables. The price to pay for this 
advantage is the introduction of a new length scale into the model: the thickness 
W of the phase-field front. For a given macroscopic problem, simulations with 
a phase-field model yield in general results that depend on the value of W. 

In the field of crystal growth, great progress towards efficient and precise 
simulations has been made by reducing this dependence on the interface thick- 
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ness [51 [3 [HI HI HD] ■ This was made possible by a detailed analysis of the model 
equations using the method of matched asymptotic expansions, which is a sys- 
tematic procedure to calculate the effective boundary conditions "seen" by the 
macroscopic transport field. Since this analysis is carried out within a pertur- 
bation approach, these boundary conditions are naturally expressed as a power 
series in W. Within the phase-field community the limit — is referred to as 
the sharp-interface limit. When the corrections due to the finite interface thick- 
ness are taken into account for choosing the model parameters, the accuracy of 
the phase-field method can be drastically improved. This procedure, which has 
been called thin-interface hmit [B], has so far been worked out only for a few 
specific physical systems. 

To be more precise, let us consider the problem of solidification, in which the 
relevant transport process that limits the growth of the crystal is the diffusion 
of heat and/or solute. Two cases are completely solved: the symmetric model, 
where the diffusion coefficients in the two phases are identical [B], and the one- 
sided model, in which no diffusion takes place within the solid phase [SI [S]. 
However, so far no method has been found to eliminate all thin-interface effects 
in the case of arbitrary diffusion coefficients in the two phases, despite some 
recent progress [71[TT1[2] (for a more detailed discussion, see [T^V 

As will be pointed out here, part of this problem arises from the fact that 
for a truly two-sided model (with different diffusivity in each phase) even the 
stationary transport problem without interface motion exhibits thin-interface 
effects. This prevents a solution of the problem by the antitrapping approach, 
which has been successful for the one-sided model [3 |9] and for the two-sided 
model with vanishing diffusion current in one phase 

In fact, such thin-interface effects are fairly general and arise in a whole class 
of problems, namely, two-phase transport through a complex structure. Exam- 
ples of relevant physical situations are the conduction of electrical current or 
heat through a two-phase material with different conductivities, the magnetic 
flux through a two-phase material with different susceptibilities, or fluid flow 
though a porous medium with variations in permeabilities. At first sight, the 
advantages of using the phase-field method in these cases are less obvious since 
the interfaces do not move and therefore the problem of front tracking does 
not arise. However, the geometrical representation of complex-shaped surfaces 
can be difficult even without interfacial motion, especially in three dimensions. 
Additionally, the use of a stationary phase-field function makes it possible to 
prescribe a given boundary condition at the interface in a straightforward man- 
ner |13[ I14j , and to impose arbitrary boundary conditions at the border of a 
physical domain of complex shape, see [151 US] and references therein. 

The problem of representing a complex interface through a diffuse boundary 
gains additional relevance because the use of tomographic methods for struc- 
ture determination becomes more and more widespread. In such methods, the 
structure representation takes the form of a matrix consisting of discrete pixels 
(or voxels in three dimensions) that contain binary or intensity data indicating 
whether a point in space is "filled" or "empty". From these data, the "true" 
structure (represented for example by discrete sharp surface elements) has to be 
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Figure 1: (Color online) Total current as function of the ratio between the 
interface thickness and the radius of the circle in the cases of direct (black 
circles), inverse (red diamonds), and tensorial (blue squares) interpolation. We 
solve the spherical inclusion problem with parameters: Mi = 1, M2 = 1/2, 
R = 0.25, h = 10'^, and p = 10^^. Lines are obtained by quadratic regression 
of the simulation data (see the main text for further details). The coefhcients 
of these regressions are contained in Tab. [1] All units are arbitrary. 

reconstructed by image analysis techniques. The phase-field method is an easy 
and robust method to obtain a smoothened representation of such data [T71 HH] . 
It could be interesting to use directly this smoothened structure for accurate 
calculation of transport processes instead of going through the additional steps 
of determining the "sharp" surface geometry. 

However, it will be shown below that thin-interface effects are also present 
in the "simple" problem of two-phase transport, even without interface motion. 
Therefore, these effects must be quantified and if possible eliminated. As we 
will demonstrate below, two effects that depend on the interface thickness are 
present: transport along the surface, and an interface resistance. In the standard 
phase-field formulation, where the transport is described by a scalar coefficient 
whose value depends on the phase field, these two effects cannot be eliminated 
simultaneously. In contrast, if the transport coefficient is allowed to become a 
tensor inside the diffuse interfaces, there are enough degrees of freedom in the 
model to eliminate both effects. 
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2 Analysis 



2.1 Problem formulation 

All the problems listed above have a common structure, namely, the flux of a 
conserved quantity is driven by a potential gradient, 

j = -M{,p)VV, (1) 

where ^ is a potential. The structure of Eq. ([T]) is standard in out-of-equilibrium 
thermodynamics: a linear relationship between a flux and a thermodynamic 
driving force (the potential gradient). The mobility coefhcient M{4>) depends 
on the phase field (/). If the transport process is electric conduction, V is the 
electrostatic potential and M the conductivity; for diffusive mass transport V 
is the chemical potential and M the atomic mobility. 

The transported quantity satisfies a conservation law, which is valid both in 
the bulk and at the surfaces, and for a time-independent solution (steady-state 
flow) reads 

V-j = 0. (2) 

The problem specification is completed by a boundary condition for the potential 
at the interfaces. We assume continuity of the potential, 

V+=V-^, (3) 

where V+ and V- are the values of the potential when the interface is approached 
from the two sides. This corresponds to a rapid exchange of the transported 
quantity between the two sides of the interface. 

Since we consider a fixed and immobile two-phase structure, the phase field 
is independent of time. We will assume that the two constant values that 
designate the two phases are </) = and (p — 1, and that the field </> varies 
between these two limits continuously through a front region of width W. For 
the sake of concreteness, the reader may have in mind a sigmoid function such 
as [1 -I- tanh(a;/W^)] /2, but the explicit form of this function is not important. 
The only hypothesis we make is that for a straight interface, the profile of (j) is 
odd with respect to the point 0=1 /2, that is, (j){x) = 1 — (j){—x) for an interface 
centered at cc = 0. This is the case in all standard phase-field models. 

2.2 Surface current 

For simplicity of exposition, it is useful to focus on a concrete example. Consider 
the conduction of electric current through a two-phase material. Then, V is the 
electrostatic potential, and M{(j)) is the phase-dependent electric conductivity. 
Furthermore, consider a straight interface normal to the x direction, centered 
at a; = 0. A potential gradient along the y direction (along the interface) is 
imposed by sandwiching the material between two parallel plates located at 
zLL/2 that are held at constant potentials ±U. Since the phase field (p and 
hence the conductivity M are constant along any line of constant x (although 
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interpolation type 


Jo 


Cl 


C2 


direct 
inverse 
tensorial 


1.752751 
1.752751 
1.752749 


0.73889X 10-^ 
-1.51884X 10-1 
-1.352X 10-3 


-0.277364 
-0.289096 
-0.358637 



Table 1: Coefficients of the second order regression J(e) = Jq + cie + C2e^ for 
the three different interpolation methods. These coefficients are obtained from 
a quadratic regression to the data indicated by the symbols of Fig. [TJ 



their values differ for different values of x), Eq. ([2]) yields a constant potential 
gradient U/L directed along the y direction. Therefore, the total current J that 
flows between the two plates is given by 

J = / dx. (4) 

Since we have considered a sample that extends to infinity, this current is clearly 
infinite. However, we will be concerned only with the excess of this current with 
respect to the sharp-interface value. The latter is obtained as the current that 
would flow if </)(a;) was a step function, that is, the space between the two plates 
is filled with material 1 of conductivity Mi for x < 0, and with material 2 of 
conductivity M2 for x > 0. This yields 

U f°° U 

J= Mi-dx+ M2- dx. (5) 
J-00 L Jq L 

The difference between these two expressions is the excess of current 5 J due to 
the variation of conductivity over a zone of finite thickness. This excess 



5 J = / [M{(l,) - Ml] -dx+ [M{<f>) - M2] - dx, 
J-00 ^ Jo ^ 



(6) 



is localized in the interface, and can therefore be interpreted as a additional 
current along the surface. It can be written as the product of the potential 
gradient and a surface conductivity 

/O poo 
[M{(j,) - Ml] dx+ [M{(l)) - M2] dx. (7) 
-00 JO 

This surface transport coefficient has two obvious properties: (i) for an interface 
profile of fixed functional form, (p{x) = f{x/W), Ms is proportional to the 
interface thickness (as can be shown by a simple change of variables), and (ii) 
it is strictly zero for any value of W if 

x; POO 

[M{4)) - Ml] dx^ [M{(j)) - M2] dx. (8) 
Jo 
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In this case, the excesses of current on the two sides of the interface exactly 
compensate. For a phase-field profile that satisfies <j>{—x) = 1 — 4'{x), this can 
be simply achieved by choosing 

M(<^) =Mi0 + Af2(l-0). (9) 
This will be called direct interpolation in the following. 



2.3 Surface resistance 

Let us now analyze again a planar interface normal to the x direction, but 
this time crossed by a steady current J± along x. In this case, the continuity 
equation immediately yields that the current is constant (independent of x). 
Then, the potential V satisfies the simple equation 



M{(j))d^V = Ji. 



Integration along x yields 



(10) 



(11) 



where V is the potential at a; = (an integration constant). In contrast, if the 
interface was sharp, the potential would simply be given by 



Voix) - V 



xJ± 



(12) 



for a; > and a; < 0, respectively. Of course, outside the diffuse interfaces, the 
slopes of V{x) are identical in Eqs. ([TT]) and ([T2|) . Therefore, the asymptote of 
the diffuse-interface expression is of the form 



v{x) - y « - 



xJ± 
Ml" 



(13) 



for X ±oo. The constants and F_ (the interface potentials "seen" from 
the region outside the diffuse interface) are readily obtained from the matching 
of this expression to Eq. pT|) . 



V+.- =V + Ji 



OO, — CX3 



1 



M{(j)) Ml, 



dx. 



(14) 



Of particular interest is the fact that these surface potentials can be different, 
in contradiction to the assumption of Eq. The difference SV = V+ — 

can be written as the product of the current J± and an interface resistance 



Rs = 



1 



1 



M{(j)) Ml 



dx 



1 



M(0) 



1 

M2 



dx. 



(15) 
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Figure 2: (Color online) Rate of convergence of the three interpolations as 
function of the ratio between the interface thickness and the radius of the disk. 
Black circles, red diamonds, and blue squares stand for direct [Eq. ([5])], inverse 
[Eq. p6|) ]. and tensorial interpolation [Eq. ((T7| ]. respectively. The physical and 
numerical parameters are the same used in Fig. [T]and Jq = 1.75275. Red and 
black solid lines are guide to eye with slope equal to 1.0, while blue dashed line 
has slope equal to 2.0. All units are arbitrary. 



This resistance is often called Kapitza resistance and has been frequently ob- 
served in experiments and simulations 19, 20, 2ll[22l[23]. Again, it is obvious 
that it is proportional to the interface thickness, and that it vanishes if the 
integral is exactly zero. This can be achieved for any value of W by the inter- 
polation 

This will be called inverse interpolation in the following. 



2.4 Tensorial mobility 

In summary, the two interface effects (surface current and surface resistance) 
can each be eliminated by a specific choice for the interpolation function of the 
mobility. Since these interpolations are mutually exclusive, it seems as if nec- 
essarily one of the two effects must remain nonzero. However, the current is a 
vector quantity, and the two effects are linked to distinct components of the cur- 
rent vector: the excess surface conductivity is relevant only for the components 
parallel to an interface, whereas the surface resistance modifies the boundary 
conditions for the normal component. Inside the two phases, where each medium 
is isotropic, the Curie principle requires to choose a scalar mobility to relate the 
current and the potential gradient. However, in the presence of an interface, 
isotropy of space is broken and a tensorial transport coefficient is permitted. In- 
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Figure 3: (Color online) Streamlines for the second geometrical configuration 
studied. The black dash-dotted, red dashed, and blue solid lines correspond to 
direct, inverse, and tcnsorial interpolation, respectively. The simulation param- 
eters are Mi = f , M2 = 0.1, R = 0.25, h = lO"^^ p = 10"^ and W = 4/i. For a 
larger version of this picture see |25| . All units arc arbitrary. 

deed, the gradient of the phase field can be readily used to define the interface 
normal n = V0/|V(/)|, which provides a second direction that is independent of 
the potential gradient. Then, we can define the transport coefficient by 

M(^) = 7\f^n (g) n + M||(l - n ® n), (17) 

where 1 is the unit tensor, with two independent interpolation functions M±{(j)) 
and M| I ((/)). If we interpolate M± according to Eq. and M|| according to 
Eq. both thin-interface effects are eliminated. Hence, for this interpolation 
the transport problem defined by Eq. ([l} becomes 

j = -M{4>) ■ VF, (18) 

with components (in two dimensions) 

jx = M,,d,V + M^ydyV, (19) 
Jy ^ M^yd^V + MyydyV. (20) 

Here, we have designated by Alij the elements of the symmetric tensor M {(()). 
The simple calculations developed above are valid only for planar interfaces. 
However, for a sufficiently smooth interface (that is, with a local radius of cur- 
vature R satisfying R 3> W), a local curvilinear coordinate system can be 
defined in which the above relations remain valid at least up to second order in 
e = W/R. It should be mentioned that a similar strategy has been used recently 
to develop efficient phase-field models for surface diffusion 
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Figure 4: (Color online) Total current as function of the ratio between the 
interface thickness and the radius of the circle for direct (black circles), inverse 
(red diamonds), and tensorial (blue squares) interpolations. The simulation 
parameters are the same of Fig. [3l All units are arbitrary. 

3 Numerical validation 

We quantify the thin-interface effects in the three different interpolations of the 
transport coefficient by solving the problem defined by Eq. ([T]) [or (fT8|) ] and Eq. 
([2]) in a simple geometrical setup. We consider a square domain T> = {{x,y) G 
[0, 1] X [0, 1]} with Dirichelet boundary conditions on the lateral edges and zero 
flux at the upper and bottom edges. More precisely, we impose V(0,y) = 1, 
V{1, y) = —1, and dyV{x, 0) = dyV{x, 1) = 0. In the center of the domain we 
place a disk of radius R. The mobility coefficient takes a value of Mi (M2) 
outside (inside) the disk, respectively. We refer to this geometrical setup as 
spherical inclusion problem. 

By combining the flux equation Eq. ([T]) [or its tensorial counterpart (fT8|) ] 
with the conservation law of our problem, i.e. Eq. we obtain an elliptic 
equation for the case of direct and inverse interpolation 

V • [M{(j)) yv] = 0, (21) 

whereas the tensorial interpolation leads to the following equation 

V ■ [M (0) • VV] = 0. (22) 

Details about the discretization and the method used to solve these equations 
can be found in the appendix. 

The thin-interface effects arising in the spherical inclusion problem are quan- 
tified by measuring the total flux at x = 1, 

I jx{l,y)dy (23) 
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and by plotting J as a function of the ratio e = W/R between the interface 
width W and the radius of the disk R, see Fig. [1] As shown in Fig. [U the three 
interpolations converge to the same value of Jq when e — ?► 0. The estimation of Jq 
is given by a quadratic regression in e of the simulation data. The coefficients 
of the regressions are listed in Table [TJ Within the truncation error [that is 
O (/i^) ^ 10~^] the three interpolations give the same value of Jq = 1.752750 ± 
10~^, as shown in the second column of Table [TJ 

In addition, we have estimated the rate of convergence of J(e) to the sharp- 
interface limit Jo for the three interpolations of the mobility. As shown in Fig. 
[21 the direct and inverse interpolations converge only linearly with e to this 
limiting value, whereas the tensorial interpolation suppress linear thin-interface 
efi^ects, which leads to a convergence that is almost quadratic in e. 

For further illustration, we consider a second geometrical configuration for- 
med by two half-disks of radius R placed in a rectangular domain with the same 
boundary condition of the spherical inclusion problem, see Fig. [3] We choose 
a small value of the ratio between M inside and outside the half-disks, e.g. 
M2/M1 = 0.1, and in Fig.[3]we plot the streamlines for this configuration. In this 
configuration, the flux is constricted in the narrow space between the two half 
disks. The difference between the three interpolations is largest for streamlines 
which are locally almost tangent to the half-disks; note the divergence of the 
different lines close to the tips of the half-disks. As shown by Eq. ([5]), choosing 
M according to the direct interpolation cancels exactly the surface conductivity 
Ms for a flat interface. Therefore, thin-interface errors affecting streamlines 
which are parallel to the boundary of variation of the transport coefficient are 
more pronounced in case of inverse interpolation (red dashed lines) than in case 
of direct interpolation (black dash-dotted lines), see Fig. [3] In Fig. [U we show 
again the values of the total current as a function of the ratio e, and again the 
tensorial interpolation performs much better than the other two, as expected. 

4 Conclusion 

We have investigated a phase-field model with a mobility tensor, in which normal 
and parallel components of the flux are interpolated with distinct functions of 
the phase fleld. Contrary to phase-field models with scalar mobilities, this 
method makes it possible to eliminate at the same time the additional surface 
diffusion and the surface resistance, which are both linked to the finite thickness 
of the interfaces. This opens the possibility to perform accurate simulations of 
two-phase transport problems with enlarged interface thickness, which can lead 
to dramatic savings in computation time. 

The complete elimination of thin-interface effects to first order in the inter- 
face thickness is also a prerequisite for the development of a quantitative crystal 
growth model for arbitrary ratio of the diffusion coefficients in the two phases. 
Indeed, a tensorial diffusion coefficient can remove both surface diffusion and 
the Kapitza resistance [12] from such models. Even if a second order asymp- 
totic analysis of a time dependent version of the tensorial problem ()17|) has 
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been recently performed [53], this methodology is still not capable of removing 
all thin-interface effects from phase-field models with a non-stationary The 
obstacle that needs to be overcome for the successful development of such a 
model is to find a coupling of the transport equation to an evolution equation 
for the phase field that yields the correct boundary conditions for the transport 
field at moving interfaces. We hope to be able to report on this problem in the 
near future. 
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Numerical Methods 

The differential operators involved in Eq. (|21l) and (|22p are approximated by 
finite differences. The domain V is discretized by an uniform square mesh with 
elements of area h"^. The fields M, V, jx, and jy are placed in different grids 
according to the standard marker and cell (MAC) method ^7\. In the primary 
grid we locate (p, M {M^x and Myy in case of tensorial interpolation), and V 
while in the two staggered grids (one for each direction) we put each component 
of the flux j. The primary nodes («, j) have coordinates 



whereas the two staggered grids are shifted of h/2 in the x (y) direction for jx 
(jy), respectively. In what follows, with the notation [. . wc mean that the 
field inside square brackets is evaluated at the position (i,j) with respect to 
the primary grid. For example, at the nodes of the x staggered grid we have 
jxij that is the value of x component of the flux between the nodes (i, j) and 
(i + 1, j) of the primary grid, i.e. [jx]i+i/2,j- 

In order to compute j with this MAC arrangement we have to specify the 
values of M at the nodes of the staggered grids. This is readily done by using 
the averaging operator on the k coordinate fi^ , that is 



X 



Hi + 1/2), 



(24) 
(25) 



y = /i(j + l/2) 



M*±i/2,, = fJ-tM^,, = (M,±ij + Af,,,) /2, 

M^,±l/2 - = (M.,J±1 + M,;j) /2. 



(26) 



(27) 



Hence, the two components of the flux j read 



jxi,j — [jx]i+l/2j- — f^x^'^ij '-^x^i,j/f^i 



(28) 



(29) 
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where h~^A~^ is the standard forward difference operator acting on the k di- 
rection. Finahy, Eq. ([2T|) is discretized by employing the backward difference 
operators /i^^A^T 

^ E ^ki^^tM^.,KV^,)^0. (30) 
k—x.y 

Evidently, the discretization of the elliptic problem arising from the tensorial 
interpolation M is more involved. In fact, two parts of Eq. (|22|) mix x and y 
derivatives, i.e. dx {MxydyV) and dy {MxydxV) . In order to guarantee second 
order accuracy and maximum compactness of the discrete stencil it is convenient 
to place the M^y component of the tensor M on a third grid, whose nodes are 
shifted by h/2 in the two directions [28]. The tensorial interpolation produces 
an additional contribution to 

[MxydyV]^^,^,^^ = fly (A4,,,,A^+A+T/,,,) /h, (31) 

and to jy 

[MxydxV\^^y^ = fi- (A4^,,,M+A+F,,,) /h. (32) 

As before, by applying backward differentiation we obtain the discrete version 
of Eq. ^ 

^{A- [^ltMxx^,J A+V,,,+fi- {Mxy^,J^Jit^t^^,3)] + 

h ^ (33) 

The MAC arrangement ensures a second order accuracy of the truncation 
error of the two elliptic problems HTJ HH]. Eq. ([5(1)) and (155)) are two linear 
systems of equations where the unknowns are the values of the potential V^.j. 
These linear systems can be easily solved through any iterative method, for 
example by the Successive Over Relaxation (SOR) [29,. To find an accurate 
solution of these problems we iterate the SOR algorithm until the maximum 
residue at the nodes of the primary grid 

p = max|[V-i], .|, (34) 
is comparable with the truncation error, i.e. p <h?. 
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